function [ BW cf1 h Ilog] = MaskingTools(volume,mask_thresh)
%% MaskingTool is a debugging function to clean and filter a contrast
% enhanced 9.4T MR image.  The steps below are used in the pipeline.
% The input volume is typically a 3D image (mandatory) around 512x512x512 in size
% with resolution 50.8 micrometers (isotropic). Binary mask, BW, is the
% primary output, but other output variables are included for debugging.
%
%   - Michael Eager,   (michael.eager@monash.edu) 
%   - (c) 2012, Monash Biomedical Imaging, Monash University, Australia

%     Copyright © 2012-2013 Michael Eager <michael.eager@monash.edu> 
%
%     This file is part of Xglom.
% 
%     This is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.


Xglompath = which('Xglom');
[Xglom_path,Xglom_filename,MRtools_ext] = fileparts(Xglompath);
addpath(fullfile(Xglom_path,'plugins'));


if nargin < 2
    mask_thresh = 2;
end



tstart=tic;
%% 1.  Use RegularizeMRI to find 2D mask of image without sample
tic;mask = RegularizeMRI(volume,499,504,1);toc

%% 2.  Apply mask to clean image
tic;cf1 = MaskVolume(volume,mask);toc

%% 3. Apply 3D Laplacian of Gaussian filter to cleaned image
% This identified gaussian blobs in the image
h = fspecial3('log',5);
Ilog = imfilter(cf1,h,'replicate');
% 31.530727 seconds.


% tic
% h = fspecial3('gaussian',10, 5);
% Inew = imfilter((Ilog.^2).*(Ilog>0),h,'replicate');
% toc
% 265.989959 seconds.


%% 4. Smooth LoG output with inhomogenityCorrection
%  Smudge gaussian blobs together
tic
corrMask=inhomogenityCorrection((Ilog.^2).*(Ilog>0),3./60);
toc
%  5.78 sec
[h x] = hist(corrMask(:),0:12/999:12);

%% 5.  Apply threshold to smoothed image and get connected components of binary image
% Binary threshold should group 'glomeruli' blobs close to each other making them the largest component. If the largest
% component is not the appropriate size include the next largest component.

CC = bwconncomp(corrMask>mask_thresh);  % mask_thresh is selected from histogram
numPixels = cellfun(@numel,CC.PixelIdxList);
[x idx]=max(numPixels);
if x < 15000000
    [x idx]= sort(numPixels);
    if sum(x(1:2)) < 15000000
        disp('Binary labels not big enough to encompass two kidney glomerular regions')
    end
    
end

%% 6. Combine the large(st) component(s) to form a binary mask
BW = ismember(labelmatrix(CC), idx);
%  9.016970 seconds.

disp('MaskingTool Elapsed time')
telapsed = toc(tstart)
%  53.97 sec


%%  Testing for clean masks of kidneys
% 

% cd ~/Monash016/RatKidney/Agilent/
% load noise_reduced_mag/cf_iso_01.mat
% load noise_reduced_mag/cf_iso_02.mat
% load noise_reduced_mag/cf_iso_03.mat
% cf(:,:,:,1)= abs(cf_iso_01);
% cf(:,:,:,2)= abs(cf_iso_02);
% cf(:,:,:,3)= abs(cf_iso_03);
% meancf = mean(cf,4);

% tic; imf = imfill(meancf.*BW2,'holes');toc
% Elapsed time is 56.070345 seconds.
% tic; imf = imfill(meancf,'holes');toc
% Elapsed time is 123.568135 seconds.
